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Abstract 

This  project  developed  the  technology  required  to  develop  electro-thermal  behavioral  models  of  microwave  circuits 
and  subsystems  Electrical  behavioral  models  themselves  are  poorly  developed  however  we  have  developed  an 
electrical  behavioral  model  that  can  be  readily  extracted  from  single  tone  measurements  or  simulations  and  capture 
baseband  effects  and  in  particular  themial  effects.  The  core  technology  is  based  on  mathematics  developed  by  the 
authors  called  the  Arithmetic  Operator  Method  which  enables  the  basic  mathematical  operations  such  as  addition, 
multiply  divide  and  subtraction  of  signals  to  be  performed  by  operating  on  their  spectra.  An  example  of  application 
of  the  method  to  the  modeling  of  a  logarithmic  amplifier  are  presented.  The  model  was  implemented  M  ATLAB 
code  and  this  code  is  available  by  contacting  m.b.steer@ieee.org.  The  code  has  been  transferred  to  two  MURIs 
(Multidsciplinary  University  Research  Initiatives)  managed  by  the  Army  Research  Office.  In  the  first  MURI  the 
work  developed  here  will  be  used  to  model  quasi-optical  power  combining  amplifiers.  In  the  second  MURI  the 
work  will  be  used  to  model  the  co-site  interference  problem. 
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1.  Introduction 

The  Arithmetic  Operator  Method  (AOM)  was  pioneered  by  Chang  and  Steer  and  has  previously  been 
described  in  two  publications,  [1]  and  [2].  It  has  been  successfully  applied  to  electronic  circuit  behavioral 
modeling  by  Gard  et  al  [3]  and  by  de  Carvalho  and  Pedro  [4],  [5].  This  project  further  developed  this 
behavioral  modeling  approach  and  implemented  it  in  MATLAB  together  with  a  worked  example  and 
comparison  to  SPICE  results. 

At  the  most  basic  level,  the  Arithmetic  Operator  Method  (AOM)  is  a  computer-aided  method  for  simulating 
the  time-domain  multiplication  of  two  signals  in  the  frequency  domain,  using  convolution  in  the  frequency 
domain.  However,  its  applications  are  not  limited  to  the  simple  multiplication  of  two  time-domain  signals 
that  one  would  find  in  a  typical  radio-frequency  mixer,  as  a  case  study  will  later  illustrate.  We  will  see  that 
AOM  can  be  used  to  model  transfer  functions  of  considerable  complexity,  as  long  as  three  prerequisites  are 
met: 

I  The  transfer  function  to  be  modeled  can  be  expressed  as  a  truncated  infinite  series  form  and  must  be  a 
mathematically  “analytic”  function  as  defined  by  Taylor  and  Mann  [6]; 

2.  The  spectral  content  of  the  input  is  known;  and 

3.  The  domain  of  the  spectral  content  (i.e.  frequencies,  but  not  amplitude  or  phase)  of  the  output  is  either 
known  or  pre-specified. 

At  first  glance,  the  last  prerequisite  might  seem  difficult  to  meet.  However,  we  will  see  that  when  the  input 
is  limited  to  a  finite  number  of  sinusoids  of  arbitrary  frequency,  then  it  is  possible  to  specify  the  domain  of 
the  spectral  content  of  the  output.  In  the  case  study,  we  demonstrate  how  AOM  can  be  used  to  model  the 
response  of  a  logarithmic  amplifier  to  a  two-tone  sinusoidal  input  signal. 


2.  Mathematical  Description 

2. 1  Mathematical  Foundation 


The  theoretical  background  underpinning  the  Arithmetic  Operator  Method  was  first  described  by 
Oppenheim  and  Schafer  [7]  in  1975.  In  chapter  10  of  their  book.  Digital  Signal  Processing,  the  authors 
describe  a  “Generalized  Superposition”  principle  that  has  extended  utility  beyond  the  usual  class  of  Linear, 
Time-Invariant  (LTI)  systems.  The  essential  properties  of  such  systems  is  that  their  operators  must  be  both 
commutative  and  associative  over  the  domains  on  which  they  operate.  Boyce  and  DiPrima  [8] 
demonstrated  that  not  all  functions  which  can  be  convolved  are  associative,  and  this  puts  some  restrictions 
on  the  form  of  the  functions.  The  functions  we  consider  will  meet  those  restrictions.  Such  systems  can  be 
modeled  using  a  linear  transformation  on  a  vector  space,  i.e.  using  matrix  algebra.  Systems  that  perform 
multiplication  in  the  time  domain,  which  corresponds  to  convolution  in  the  frequency  domain,  exhibit  the 
properties  required  for  “generalized  superposition”  -  even  though  they  are  inherently  non-linear.  The 
resulting  linear  transformation  implements  the  essence  of  the  Arithmetic  Operator  Method.  We  present 
here  the  development  of  the  transformation. 

2.2  Time  and  Frequency  Domain  Descriptions  of  Signals  of  Interest 


Let’s  begin  with  two  time-domain  signals,  x(/)  and  z{t)  ,  both  of  which  are  composed  of  a  finite  number 
of  sinusoids  of  arbitrary  frequency  content.  For  the  purposes  of  the  development  here,  we  will  assume  that 
X  and  z  have  spectral  content  at  the  same  frequencies,  although  with  different  amplitudes  and  phases.  This 
is  not  restrictive,  since  if  x  and  z  were  composed  of  sinusoids  of  different  frequencies,  we  could  describe 
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both  using  a  common  (i.e.  union)  set  of  frequencies  with  some  amplitude  coefficients  in  x  and  z  both  set  to 
0,  Thus,  the  development  that  follows  here  for  xit)  will  also  apply  to  z{f) .  We  can  express  x{t)  as  e 
sum  of  a  set  of  sinusoidal  signals  as  follows: 


x(t)  =  Xa,-„(0  =J,\x„\cos(.a)j+<p„)  =^\x„ 


(1,1) 


Here,  we  have  expressed  each  sinusoid  using  the  cosine  function  with  frequency  Ct)„  -  27tf„  ,  amplitude 

lx  I  and  phase  0  .  Following  that,  we  have  used  Euler’s  equation  =  COS0 +7  Sin0)  to  express 

the  cosine  function  as  the  sum  of  two  complex  exponentials.'  We  can  further  refine  this  result  as  follows: 


-;(W„  <+<?„) 


«=0 

N 


n=0 

N 


H=0  »=0 


(1.2) 


Here  we  have  again  used  Euler’s  equation  to  express  the  amplitude  and  phase  as  a  complex  number  pair, 
X  +  iX  ■  If  we  now  take  the  Fourier  transform  of  x{t}  using  the  final  form  in  equation  (1.2)  above, 

—  J  HI  ■ 

we  will  have: 


XiCO)  =  F  {x(0}  =  F  hY,[X„r+jX„iy‘‘"'  +l^[Xnr  " 
I  »-o  »=o 


(1.3) 


/i=0 


Now,  from  the  definition  of  the  inverse  Fourier  transform,  we  have: 


‘  This  form  will  allow  us  to  use  a  phasor  form  for  the  cosine  function  later  in  the  development. 
However,  it  must  be  stressed  that  since  we  are  developing  a  model  for  a  non-linear  system,  we  cannot 
make  the  usual  linear  phasor  transformation  x(0  =  Rc[A^(fi))] ,  and  then  use  X{03) .  Our  representation 
is  exact,  although  we  will  get  the  benefit  of  phasor  arithmetic  through  our  formulation. 
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x(0  =  F-’{X(fi))} 

+  ;Xj[2;r5(ffl-®j]+ii[X„.  -7Zj[2;r5(®+0j] 

”  i  [  «=0  »=»  •' 

-1,1  »=0  ”="  J 

=T{iiK + 

(1.4) 

Thus  we  see  that  there  is  a  factor  of  2JI  that  appears  when  we  follow  the  strict  definition  of  the  Fourier 
transform.  This  same  factor  is  then  eliminated  in  the  inverse  Fourier  transform.  For  the  purposes  of  our 
development  we  will  ignore  this  domain  conversion  factor  as  we  go  forward,  since  the  ultimate  result  of 
our  efforts  will  be  a  time  domain  signal.  With  this  simplification,  we  observe  that  the  frequency  domain 
constants  for  each  5{(0-(0j  component  are  simply  |  (X„,  +  jX,,-,  )  =  jX„  for  components  on  the 
right  hand  side  of  the  spectrum.  On  the  left  hand  side  of  the  spectrum,  each  5(0  +  ®„)  will  have  a 
multiplicative  constant  of  i(X„.  -  jX„,)  =  ^(X„,  +jXj  =  ,  where  X’,,  denotes  the  complex 

conjugate  of  X„ .  Of  special  note  is  what  happens  for  the  case  of  DC,  when  n  =  0  .  By  definition,  DC 
corresponds  to  the  origin  in  the  frequency  domain,  so  there  is  no  imaginary  part  to  the  constant  X^ ,  i.e. 
X^  =  Xor .  This  can  also  be  seen  by  going  back  to  equation  (1.2),  substituting  0  for  n,  and  working 
through  the  complex  arithmetic.  Note  also  that  there  is  no  scale  factor  of  ■j  on  the  DC  temi.  We  can  now 
rewrite  X{CO)  from  equation  (1.3),  dropping  the  factor  of  2k  and  taking  account  of  the  coefficient  scale 
factors: 


To  recap,  the  first  summation  is  a  set  of  frequency  domain  components  on  the  right  hand  side  of  the 
spectrum,  the  second  is  a  set  of  frequency  domain  components  on  the  left  hand  side  of  the  spectrum,  and 
we  carry  the  DC  term  separately.  Now,  since  we  will  be  dealing  with  convolution  in  the  frequency  domain 
shortly,  it  will  be  convenient  to  designate  the  components  in  the  left  hand  side  of  the  spectrum  with 
negative  indices,  so  we  have: 


X(o})  =  iYx„5(co-co„)  +  X,8ico)+^f^x:„5ia>+o}.„)  (>  6) 

/i=l 

. .  .or,  in  an  even  more  compact  form: 

X{co)  =  X„5  (ffl)  +  f  (to + £0.„ )}  (1.7) 

n=l 
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Now,  we  can  express  Z((0)  similarly,  i.e.: 


Z(co)  =  Z„S(co)+i'Z(Z,S(c^-co„)  +  z:„5(co+o)_J} 

n=l 


(1.8) 


2.3  Convolution  with  Reduction  of  the  Frequency  Domain  to  Positive  Frequencies 


At  this  point  we  can  consider  the  effect  of  defining  the  time  domain  operation  y(0  -  x(t)  ■  z{l) .  Since 
time  domain  multiplication  corresponds  to  frequency  domain  convolution,  we  have 
y{t)  =  X{t)  ■  z{t)  <^Y{ai)  =  X (ffl)  *  Z(ffl)  ,  so 

r(<B)  =  |Xo5 +  * 

L  ■'  (1,9) 

|z„5  (  ffl)  +  i  X  (®  -  ®»)  + 

This  is  a  cumbersome  equation,  but  the  effect  of  convolving  X{(0)  with  Z{a>)  will  be  a  set  of 
8  functions  with  frequencies  ranging  from  -2(0 (at  summation  index  -2N)  to  2(0^  (at  summation  index 
2N).  We  will  thus  have: 


2N 

Y{(0)  =  ZoS  (co)  +  i  £  {Z„5(£d  -  fo„ )  +  +  ffl_„ )} 


/inl 


(1.10) 


If  we  examine  the  form  of  Y{(0) ,  we  see  that  the  left  hand  side  of  the  spectrum  has  the  same  amplitude 
and  opposite  phase  of  the  right  hand  side  of  the  spectrum.  It  thus  conveys  no  information  above  and 
beyond  that  available  from  the  right  hand  side  of  the  spectrum  alone,  and  so  we  can  reduce  the 
computational  complexity  involved  in  carrying  out  frequency  domain  convolution  by  eliminating  the  left 
hand  side  of  the  spectrum  in  the  convolution  operation.  A  consequence  of  doing  this  is  that  we  must 
recover  the  amplitude  and  phase  of  each  component  of_)'(f)  from  the  coefficients  of  the  right  side  of  the 

spectrum,  i.e.  given  that  _v(0  has  the  form: 


2N 

y{t)  =  Xlz»  1  cos(£»„f +t)+yo 

H=I 


(1.11) 


Then  without  adjustments,  we  would  have,  for  each  n. 


and: 


0„=tan-' 


(1.12) 


(1.13) 
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The  phase  is  unaffected  by  our  dropping  of  half  the  spectrum,  but  in  order  to  recover  the  correct  amplitude, 
equahon  (I  12)  must  be  scaled  by  a  factor  of  2.  In  the  matrix  formulation  that  we  are  developing,  we  take 
account  of  this  factor  of  two  in  the  formulation  of  the  matrix,  leaving  the  multiplying  vector  unsealed. 

Since  the  number  of  convolution  operations  is  reduced  geometrically^,  while  the  additional  operations 
increase  linearly,  this  is  a  good  tradeoff.  As  we  proceed  with  the  development  of  the  method,  we  will  show 
frequency  domain  convolution  only  for  the  right  hand  side  of  the  spectrum  with  the  expectation  that  the 
amplitude  and  phase  of  y{t)  are  determined  as  above.  With  the  spectrum  limited  to  the  right  hand  side,  we 

define  the  spectrum  forx{t)  from  equation  (1.1)  as: 


Xim)  =  X„5(co) + i  ) 


(1.14) 


Nowz(0  would  be  equivalent  except  for  the  fact  that  we  wish  to  account  for  the  output  scale  factor  in  its 
formulation,  so  we  will  multiply  all  frequency  components  by  2  to  take  this  into  account: 

Z(o)  =  Z,5(o)+2C||;z„5(ffl-©„) 

(1.15) 

=  Zo5(ffl)+5^Z„5(ffl-ffl„) 

«=1 


Now  let’s  consider =  x(0  -2(0  given  this  redefinition  in  the  frequency  domain.  We  have 
y{t)  =  x(t)  ■  z{t)  r((B)  =  A'(ffl)  *  Z(,a)) ,  so: 

Y{a)  =  Xico)*Z{co)=  £  Xico,)Zico-(o,)(\.\6) 

Alla)^ 

Expanding  this  by  substituting  equations  (1 . 1 4)  and  (1 . 1 5)  into  equation  (1 . 1 6),  we  have: 


^  Consider,  for  example,  a  case  where  x  and  z  arc  both  composed  of  tones  of  10  different  frequencies. 
This  corresponds  to  21  frequency  domain  locations  (after  including  DC)  that  must  be  convolved  with  21 
other  frequencies,  or  computational  complexity  on  the  order  of  21^  =  441 .  Restricting  the  frequency 
domain  to  the  right  hand  side,  we  need  convolve  only  1 1  frequency  domain  components  for  a 
computational  complexity  on  the  order  of  1 1^  or  121.  The  ratio  is  just  under  4:1  or  2^:1,  indicating  the 
geometric  reduction  in  computation. 
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Y{o})  =  Xio})*Zico)=  X  X(ffl,.)Z(ffl-£B,)  = 

JZo5(®-®J+Z,5(®-®*-®,)+Zj5(ffl-®*-®2)  +  ... 


AllcOt 


(1.17) 


{Xo5(®*)+iZ,5(®* -(0,)+^X,5i(0^ -o),)+...^X^5i(0,-co^)]- 
[Zo5(ffl-®i)+Z,5[®-(®j  +®,)]+Z25[©-(®*  +ffl2)]+-l 


+Z/,5[a)-(a),  +  a)A0] 


In  the  equations  above,  we  make  only  one  stipulation  -  that  6)q  —  0 .  We  can  now  make  some 
observations  on  the  frequency  domain  convolution  above.  First,  for  a  total  of  N  distinct  input  tones  or 
frequencies  -  all  in  the  right  hand  side  of  the  spectrum  -  the  output  Y(CO)  has  components  at  DC  and  2N 
positive  frequencies,  the  maximum  radian  frequency  being  26)^  .  Second,  exhibiting  the  usual  sifting 
property  of  the  <5  function,  for  each  distinct  value  of  the  frequency  index  k,  only  one  term  in  X(COf.)  will 
be  non-zero.  This  term  will  then  form  a  set  of  products  with  all  of  the  components  of  Z  (Ct)  ~  ft)  J ,  which 
retains  the  independent  variable  ft) .  Put  another  way,  each  sifted  term  in  X {(Of. )  forms  product  terms 
with  all  of  the  terms  in  Z{CO)  ,  but  shifted  in  the  frequency  domain  to  the  frequency  dictated  by  the 
argument  to  the  5  functions  for  Z(ft))  called  for  in  equation  (1.17).  In  the  course  of  the  convolution 
operation,  there  may  be  a  number  of  non-zero  products  for  different  values  of  the  frequency  index  k  for 
which  the  5  functions  inZ(ft))  happen  to  evaluate  to  the  same  frequency.  The  convolution  operation  calls 
for  these  terms  to  add  linearly,  conforming  to  the  notion  of  generalized  superposition. 

To  illustrate  this,  consider  the  cases  of  k  =  \  and  k~2. 


When  A:  =  1 ,  we  have: 


Zo5(®— ffl,)  +  Z,5[®— (®|  +ffl|)]  + 


Z^5[a-{(o^  +®2)]+'-+ 

Zo(5(ffl  -  ®| )  +  Z|5(®  -  2®| )  + 

Z^5{(£>-{(o^  +(o^)]+...Z^8[0}-{o)^  +®;^)]) 


(1.18) 


And  when  k  ~2  ,we  have: 
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(1.19) 


Z|,5  (tt)  -  (Dj) + Z,5[g)  -  (fiij  +  ®i )] + 


■  Zj5[(»-(fi92  +  02)]+-  + 
Z„5[ct)-(ct)2 +  (»;,)] 


^  ^  |z25[ffl -2^2  )]+•••  ■^a-5[£o  -  (ffl, +0)^)])} 


Consider  the  frequency  component  (0^+(0^,  which  is  indicated  by  any  terms  containing 

-  (ffl|  +  ©2  )]  in  equations  ( 1 . 1 8)  and  ( 1 . 1 9).  From  equation  (1.18),  there  is  a  term  at  CD,  +  CDj  with 
the  coefficient  \X^Z^,  while  from  equation  (1.19),  there  is  a  term  at  ta,  +  COj  with  the  coefficient 
^  XjZ, .  In  general,  the  terms  in  each  Y (CO)  |  j  will  be  composed  of  spectral  components  at  different 

frequencies  -  each  coefficient  having  a  different  argument  to  its  S  funetion.  This  is  not  the  most  useful 
form  to  implement  using  computed-aided  methods,  as  we  will  now  see. 

2.4  Formulation  of  a  Matrix  Method  for  the  Frequency  Domain  Convolution 


Let’s  examine  the  result  of  performing  the  complex  multiplication  required  to  create  one  of  the  components 
of  .  Continuing  with  terms  that  appear  at  the  frequency  (0^  +  in  7 (fO) ,  let’s  examine  what 

happens  when  we  perform  the  complex  multiplication  of  X,  and  : 

^1^2  ”  (-^Ir  ■*"  )  '  (^2r  J^2i)  ~  ^\r^2r  “  ^\i^2i  ^  ^ 


The  result  we  obtain  is  a  complex  number  composed  of  four  product  terms  from  the  complex  number  pairs 
for  Xj  and  .  Of  particular  note  here  is  that  there  are  real  and  imaginary  parts  of  the  product  that  will 
combine  with  other  product  terms  occurring  at  ft),  +  (d^  to  create  real  and  complex  parts  of  7 (O))  .  Hence 
this  complex  multiplication  can  be  viewed  as  forming  two  real-valued  linear  equations.  If  we  do  this,  we 
can  formulate  the  product  in  equation  (1 .20)  in  the  following  matrix  form: 


Xl, 

^2/ 

.^ll 

(1.21) 


Here  we  use  the  primes  on  Y  merely  to  indicate  that  this  product  forms  only  a  part  of  the  output  7 (ti))  at 
the  frequency  In  order  to  get  the  complete  response  of  7 (6))  at  the  frequency  CO^+CO^,  it  will 

be  necessary  to  evaluate  each  7(fi))|jt  at  the  frequency  Ct),  +6)2  and  form  the  sum  as  suggested  by  the 
convolution  summation.  In  other  words, 

]r(dD,+602)  =  r(CD,+CD2)|*„  +  y(dD,  +  ffl2)h.2+-+J'(®l+®2)|*=2V  ('-22) 
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Equation  (1.22)  suggests  extending  the  matrix  multiplication  form  introduced  in  equation  (1 .21)  so  that  a 
simple  matrix  multiplication  will  yield  a  sum  of  products  that  forms  the  complete  response  for  Y{(0)  at  a 

particular  frequency.  As  an  illustrative  example,  consider  for  a  moment  that  Y(0)^  +(^2)  is  limited  to  just 

terms  in  A:  =  1  and  A:  =  2  .  Referring  to  the  product  terms  in  equations  (1 . 1 8)  and  (1 . 1 9),  then,  the  matrix 
form  of  equation  (1.21)  would  be  extended  to  add  extra  columns  to  the  matrix  and  extra  rows  to  the  vector: 


y; 

[(CD,+CO,) 

Y 


2., 

iXr 

^2. 

.22,- 

(1.23) 


This  formulation  can  be  generalized  by  adding  two  rows  to  the  matrix  for  each  unique  output  frequency  in 
7(C0)  ,  and  it  can  be  generalized  for  each  additional  input  frequency  by  adding  two  columns  to  the  matrix 
and  two  rows  to  the  vector  element.  The  matrix  thus  formulated  has  been  termed  the  Spectrum  Transform 
Matrix  by  Chang  and  Steer  [2].  In  matrix  algebra  notation,  the  spectrum  transform  matrix  is  denoted  , 
the  spectral  vector  of  frequency  components  for  z(^)  is  denoted  as  Z  ,  and  the  basic  operation 
y{t)  =  x(0  •  z(t)  can  now  be  written  with  its  frequency  domain  counterpart  as: 

y(t)  =  x(t)  •  z(i)  o  Yico)  -  XiCO)  *  Z(fi))  =  T,Z  (1 .24) 


2.5  Generalization  of  the  Spectrum  Transform  Matrix  to  Arbitrary  Transfer  Functions 

The  Spectrum  Transform  Matrix  as  described  in  the  previous  section  implements  the  frequency  domain 
convolution  for  the  basic  time  domain  multiplication  operation  y(t)  =x(t)'  z(t) .  However,  it  is 
possible  to  further  the  utility  of  the  Arithmetic  Operator  Method,  as  we  will  now  describe.  Consider  y(t) 

as  an  arbitrary  polynomial  transfer  function  of  x(/)  ,  for  example 

M 

=  ^0  +  a^X  +  a^X^  + ...  +  a^^X^  =  closely,  we  can  write 

m=0 

this  is  as  y{t)\„,^2  ~  ^2^^  “  a^X'X.fxom  which  it  follows  that  a2X'X  .  We  can  then 

substitute  the  matrix  product  T^X  into  the  term  for  m=3  and  obtain  a^X  •  X^  4-^  (Tt . 
In  general, 

=  =  (1.25) 

rtjs:0  m=0 

With  this  formulation,  a  word  is  in  order  about  a  convention  that  has  been  adopted  into  the  algorithm  for 
creating  the  spectrum  transform  matrix  as  described  by  Chang  and  Steer  [2].  Referring  back  to  equations 
(1.12)  through  (1 .15),  we  see  that  we  accounted  for  a  lost  amplitude  scale  factor  of^by  multiplying  the 

non-zero  components  of  Z{CO)  by  2,  leaving  a  scale  factor  of  j  in  the  Spectrum  Transform  Matrix.  Since 
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x{t)  and  z{t)  were  both  composed  of  sinusoids,  we  chose  arbitrarily  to  scale  the  spectral  vector  up  by  2. 
However,  as  we  generalize  the  formulation  of  the  Spectrum  Transform  Matrix,  the  spectral  vector  will 
generally  remain  a  set  of  sinusoids,  but  the  Spectrum  Transform  Matrix  will  not.  In  the  interest  of  keeping 

the  Spectral  Vector  description  as  simple  as  possible,  we  will  move  the  y  scale  factor  from  the  Spectral 
Vector  to  all  of  the  non-DC  terms  in  the  Spectrum  Transform  Matrix  .  By  convention,  then,  the 

Spectral  Vector  X  will  always  be  composed  of  the  unsealed  phasor  coefficients  for  each  sinusoidal 
component  of  the  input.  This  scale  factor  is  termed  £  in  the  paper  by  Chang  and  Steer  [2],  which  also 
gives  a  detailed  algorithm  for  construction  of  the  Spectrum  Transform  Matrix. 

2.6  Streamlined  Method  for  Developing  Entries  in  the  Spectrum  Transform  Matrix 


Now  that  we  have  seen  why  the  Spectrum  Transform  Matrix  construction  can  be  useful,  we  need  to  find  a 
streamlined  way  to  construct  it  -  one  that  is  simpler  than  writing  out  all  of  the  convolution  products  and 
evaluating  terms  in  F((W)  as  described  in  equation  (1.22).  We  will  see  in  a  moment  that  we  can  use  a 
simple  algorithm  to  tabulate  all  of  the  entries  in  the  Spectrum  Transform  Matrix.  This  tabular  form  of  the 
entries  in  the  Spectrum  Transform  Matrix  is  termed  the  Spectrum  Mapping  Table  by  Chang  and  Steer  [2]. 
It  is  derived  by  anticipating  all  combinations  of  the  input  tones  (including  intermodulation  distortion 
products)  that  will  result  in  incremental  phasor  components  in  the  set  of  output  frequencies  that  we  have 
anticipated  in  advance.  (Recall  that  anticipating  the  domain  of  the  spectral  content  of  the  output  is  one  of 
the  prerequisites  to  using  AOM.)  To  do  this,  we  use  a  mathematical  tool  called  the  Basic  Intermodulation 
Product  Description  (BIPD)  table. 

2.6.1  Basic  Intermodulation  Product  Description  Table 


The  Basic  Intermodulation  Product  Description  (BIPD)  table  is  a  list  of  frequencies  which  defines  the 
anticipated  frequency  domain  of  the  output  signal  y(t)  .  It  is  created  by  determining  a  set  of  frequencies 
of  interest  given  the  form  of  the  input  signal.  Often  the  input  signal  is  a  non-commensurate  two-tone 
signal,  so  we  will  use  that  for  the  purposes  of  discussion  and  illustration.  However,  BIPD  tables  can  be 
developed  for  input  signals  of  an  arbitrary  number  of  non-commensurate  signals.  As  an  example,  consider 
the  case  of  a  down-conversion  mixer  where  a  narrowband  RF  signal  and  a  local  oscillator  signal  are  to  be 
mixed  through  an  intentional  non-linearity  that  is  known  to  have  non-linear  behavior  up  to  the  third  order. 

More  formally,  the  input  is  a  signal  of  the  form  x(t)  =  AT,  COSCO^t  +  X^  COSCO^t ,  where  we  are  assuming 

2  3 

that  CO^  >  ,  and  the  mixer  has  a  non-linear  behavior  that  can  modeled  as  y-aQ+a^x-\-  a^X  +  a^x  . 

Then  because  of  the  third-order  non-linearity,  we  can  expect  to  see  harmonic  distortion  products  at  Soijand 
3G)2  as  well  as  intermodulation  distortion  products  at  2ft}j  —  (O^  and  —CO^.  In  a  down-conversion 
mixer,  we  are  not  generally  concerned  with  third-order  harmonics,  since  they  can  be  filtered  easily  through 
a  bandpass  filter.  However,  the  intermodulation  distortion  products  just  mentioned  are  of  great  concern 
because  they  surround  the  two  fundamental  frequencies  of  interest,  ft)j  and  (Oj ,  and  often  fall  within  the 
passband  of  filters  for  the  output  y.  Since  they  are  of  great  concern  to  us,  we  will  wish  to  use  AOM  to 
identify  their  possible  effect  on  the  output  signal.  The  second  order  non-linearity  will  yield  second  order 

haraionics  at  2g)j  and  ICO^  along  with  intermodulation  products  at  CO^  -  6),  (the  desired  down-conversion 
result)  and  +  tD, .  We  have  thus  identified  9  output  frequencies  in  this  example  in  which  we  are 
interested:  0  (DC),  ft)j ,  ^2  +  »  ^^2  “ 

Basic  Intermodulation  Product  Description  frequencies  of  interest  to  us.  Thus,  the  Spectrum  Transform 
Matrix  we  create  for  the  problem  should  have  a  set  of  two  rows  and  columns  corresponding  to  each 
frequency.  As  noted  in  the  paper  by  Chang  and  Steer,  since  the  DC  term  of  the  output  will  have  no 
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imaginary  part,  we  simply  delete  the  second  row  and  column  from  the  matrix  after  executing  the  algorithm 
to  create  it.  Thus  for  the  case  of  DC  and  8  non-zero  frequencies,  the  Spectrum  Transform  Matrix  will  be 
17x17. 

We  can  describe  the  BIPD  frequencies  by  means  of  a  table  that  describes  the  relationships  among  the 
frequencies  using  a  vector  algebra  form.  If  /  and  are  the  two  non -commensurate  frequencies  that 

comprise  the  input  signal  x(t) ,  i.e.  (0^  =  and  CO^  =  ,  then  we  can  describe  the  BIPD 

frequencies  by  means  of  the  following  algebraic  weighting  table: 

BIPD  Table 


kx _ nj _ n2 

0  0  0 

1  1  0 

2  0  1 

3  -1  1 

4  1  1 

5  2-1 

6-12 

7  2  0 

8  0  2 

Here,  is  a  mnemonic  index  that  we  use  to  distinguish  different  frequencies.  By  convention,  we  denote 
DC  by  =  0  and  assign  the  non-commensurate  tones  to  the  values  of  corresponding  to  the  description 
we  gave  for  x(f)  above.  With  this  construction,  we  can  now  assign  the  non-commensurate  input  tones  a 
unit  vector  algebraic  description  in  terms  of  the  frequency  weighting  factors  and  .  Doing  so  allows 
us  to  describe  the  harmonic  and  intermodulation  products  as  weighted  sums  of  the  non-commensurate  input 
frequencies.  Note,  however,  that  beyond  DC  and  the  non-commensurate  input  frequencies,  there  is  no 
functional  relationship  between  A^^and  the  harmonic  or  intermodulation  tone  it  describes.  The  index  is 

simply  a  tracking  mechanism  for  these  tones.  With  this  table  in  hand,  we  now  know  the  domain  of  the 
spectrum  of  the  output  -  one  of  the  prerequisites  at  the  outset.  That  is,  we  now  know  the  size  of  the 
Spectrum  Transform  Matrix  .  However,  we  have  not  yet  determined  the  range  of  the  output,  which  is  to 
say,  the  contents  of  the  Spectrum  Transform  Matrix,  This  comes  from  the  Spectrum  Mapping  Table,  to 
which  we  now  turn  our  attention. 

2.6.2  Spectrum  Mapping  Table 


The  rationale  for  describing  the  BIPD  frequencies  in  vector  form  becomes  apparent  when  we  confront  the 
prospect  of  trying  to  determine  the  output  at  each  of  the  frequencies  in  the  table.  Consider  again  the 

polynomial  transfer  function  of  the  output  for  our  example:  y  —  Cl2X^  +  .  Since  in  the 

frequency  domain  this  polynomial  corresponds  to  repeated  convolution  operations,  we  must  allow  for  the 
possibility  that  any  of  the  frequencies  in  the  BIPD  table  -  including  the  harmonic  and  intermodulation 
products  -  will  be  further  mixed  within  the  mixing  apparatus  to  produce  additional  phasor  components  at 
each  BIPD  frequency.  We  must  specifically  be  alert  to  look  for  instances  of  any  entry  in  the  BIPD  table 
being  subtracted  from  any  other  entry  where  the  resulting  frequency  maps  back  onto  one  of  the  BIPD 
frequencies.  We  must  determine  every  instance  where  this  does  occur,  because  the  output  y  at  each 
frequency  in  the  BIPD  table  will  be  composed  of  a  sum  of  phasor  products  of  exactly  this  sort.  For 
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example,  consider  the  result  of  subtracting  the  tone  at  frequency  index  =  3  (produced  from  the  second- 
order  mixing)  from  the  tone  at  frequency  index  k^=6,  (produced  from  third-order  intermodulation): 

Spectrum  Mapping 


kx _ Qi _ 02 

-3  1  >1 

_6 _ J ^ 

Sum  0  1 


The  result  of  this  additional  frequency  convolution  is  a  phasor  product  term  that  maps  to  the  same  output 
frequency  as  the  non-commensurate  fundamental  tone  .  There  may  be  several  phasor  component 
products  generated  in  this  fashion  that  map  onto  the  output  BIPD  frequencies.  The  enumeration  of  these 
frequencies  forms  the  Spectrum  Mapping  Table.  Unlike  the  simple  illustrative  example  here,  the  Spectrum 
Mapping  Table  is  constructed  systematically  by  fixing  the  output  BIPD  frequency  and  then  finding  all 
vector  combinations  of  the  input  BIPD  frequencies  -  including  vector  subtraction  combinations  -  that  map 
to  the  output  BIPD  frequency.  Each  entry  in  the  table  will  contain  5  pieces  of  information:  The  frequency 
index  k ,  of  the  output  BIPD  to  which  the  phasor  product  is  being  mapped,  the  frequency  indices  of  the  two 
input  BIPDs  being  mixed,  and  the  signs  indicating  frequency  addition  or  subtraction.  For  example,  if  we 
denote  our  two  input  signals  using  frequency  indices  k^  and  k^  and  denote  the  signs  of  the  vector  operation 

as  and  ,  then  the  entry  in  the  Spectrum  Mapping  Table  for  the  previous  example  would  be: 

Spectrum  Mapping  Table  Entry 

ky _ foe _ sx _ kz _ sz 

2  3-161 

It  must  be  emphasized  that  the  Spectrum  Mapping  Table  can,  and  often  does,  contain  more  than  one  entry 
with  the  same  value  for  k^  ,  this  being  the  result  of  many  different  forms  of  intermodulation  taking  place 

among  the  BIPD  frequencies.  If  we  now  refer  back  to  equation  (1.20),  we  can  formulate  the  phasor  result 
of  the  operation  specified  in  the  Spectrum  Mapping  Table  Entry  in  terms  of  the  BIPD  nomenclature: 

K  P.Z^j)= 


In  the  spirit  of  equations  (1.21)  through  (1.23),  we  can  restate  the  basic  formulation  of  the  Spectrum 
Transform  Matrix  entries  using  the  result  in  equation  (1.26).  We  will  have: 


It  is  important  to  note  the  primes  on  Y  in  equations  (1.26)  and  (1 .27)  and  emphasize  that  we  have  shown 
the  result  of  the  inclusion  of  only  one  BIPD  mixing  product  here.  As  the  paper  by  Chang  and  Steer  [2] 
notes,  the  construction  of  the  Spectrum  Transform  Matrix  is  done  by  first  initializing  all  entries  to  zero  and 
then  adding  to  the  matrix  based  upon  the  specifications  of  the  Spectrum  Mapping  Table.  Once  all  of  the 
entries  in  the  Spectrum  Mapping  Table  have  been  read  and  updates  made  to  the  Spectrum  Transform 

Matrix,  the  Spectrum  Transform  Matrix  is  now  complete.  To  begin  using  T^to  perform  frequency 


CO 

1 _ 

Zk.r 

Aj. 
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domain  convolution  operations,  all  that  remains  is  to  construct  the  spectral  vector,  Z  or  X  (which  are  the 
same  in  the  example  in  this  section,  although  we  used  the  Z  variable  for  clarity  in  explanations).  Z  must 
have  the  same  number  of  rows  as  T^,  which  in  our  example  is  17.  The  construction  is  straightforward: 

Given  that  x(0  =  COS®/  +  COSCO^t ,  Z  is  a  vector  containing  a  zero  entry  (for  DC)  followed  by 
the  phasor  values  for  its  two  non-commensurate  tones,  followed  by  a  set  of  padding  zeros  corresponding  to 
the  harmonic  and  intermodulation  terms.  Thus  we  have: 


0,0, 0,0, 0,0, 0,0, 0,0,0, of 


(1.28) 


Once  the  series  of  matrix  multiplications  corresponding  to  each  polynomial  term  has  been  carried  out  as 
called  for  by  the  form  for  y  in  (1.25),  the  result  is  a  vector  of  phasor  coefficients  for  each  polynomial  term 
in  y.  Thus  the  time  domain  y{t)  can  be  recovered  by  simply  converting  the  phasor  coefficients  for  each 
polynomial  product  into  amplitude  and  phase  terms  and  forming  the  cosine  function  with  frequency 
Infi,  as  determined  by  the  BIPD  table,  and  then  summing  the  cosines  that  correspond  to  each  polynomial 

product. 


3.  AOM  Case  Study -A  Logarithmic  Amplifier 

Consider  the  logarithmic  amplifier  shown  in  Figure  1.  Such  circuits  are  most  often  used  in  instrumentation 
systems  because  of  the  wide  range  (many  orders  of  magnitude)  of  inputs  to  which  they  can  respond  (see, 
for  example,  [9]  and  [10]).  Its  transfer  function  can  be  written  as  an  infinite  series  expansion  relating  the 
output  to  the  input  x(0 .  The  Arithmetic  Operator  method  is  used  to  solve  the  expansion  in  the 

frequency  domain. 


The  following  assumptions  are  made  about  the  circuit  elements: 

1 .  Ideal  Bipolar  Transistor 

2.  Ideal  op-amp:  The  voltage  between  terminals  tends  to  zero,  and  the  current  into  the  op-amp  is  zero.  Vo 
=  -Vbc,  and  Id  =  Ii. 


From  the  transistor  constitutive  relations: 
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=  h  [exp(^)  - 1]  □  exp(^)  (1 .29) 

where  is  the  reverse  saturation  current  of  the  transistor,  and  V-j-  is  the  thermal  voltage.  Rearranging 
terms, 


^  =  exp(^) 

(1.30) 

ln(4f)  =  ^ 

(1.31) 

(1.32) 

Now,  /.  =  4,  =.^,so 

V,=-VrH^) 

(1.33) 

V.=-Vr\n{-k) 

(1.34) 

For  use  with  the  Arithmetic  Operator  Method,  this  transfer  characteristic  should  be  represented  in  terms  of 
the  basic  mathematical  operators  -  addition,  subtraction,  multiplication,  and  division.  Three  common  forms 
exist  for  expanding  the  natural  logarithm  [11]: 

ln(x)  =  (x-l)-i-(x-l)'+i(x-l)'+...  (0<x<2)  (1.35) 

ln(x+l)  =  x-jx^+jx’-:[x‘'+...  (-l<x<l)  (1.36) 

ln(x)  =  2[fji  +  i(^)^+i(fj[)^+i(^)’  +  ...  (x>0)  (1.37) 

For  the  purposes  of  this  paper,  the  circuit  equation  was  manipulated  to  fit  the  form  of  (1.36).  The  input  can 
be  written  as 

+  -where  |^|<1  (1.38) 

Substituting  (1.38)  into  (1.34)  results  in 

^  and  this  can  be  expanded  as:  (1 .39) 

K  =  -^T  [in  (x* )  +  InCi;)  +  ln(l  +  ^)]  (1 .40) 

Of  the  terms  in  (1.40)  the  first  two  are  constant  offsets.  Only  the  last  term  has  frequency-dependent 
characteristics.  This  last  term  can  be  represented  by  its  linear  series  expansion  shown  in  (1 .36)  and  solved 
with  the  Arithmetic  Operator  Method. 

The  comparison  was  performed  for  an  input  signal  consisting  of  the  sum  of  two  sinusoids: 

V.  =  Aq+Ai  sin(2;r/j^)  +  A2  sin(27r  fjt)  ,  which  in  the  form  of  (1.38)  is: 

Vi=  A^[\  +  ^s:m{2Ttffy  +  ^sm{2nf^t)]  (1.41) 

A  A 

The  sine  form  was  chosen  in  order  to  assure  convergence  at  the  beginning  of  a  transient  simulation  in 
Spice.  Since  sine  and  cosine  are  shifted  90  degrees  in  phase,  that  is,  sin  (Ot  =  COS,{COt  —  90°) ,  we  take  the 
phase  into  account  by  designating  the  amplitudes  of  the  non-commensurate  input  tones  above  in  the 
spectral  vector  as  pure  imaginary  terms.  Parameter  values  are  shown  in  the  Table  1 : 
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Circuit  and  Input  Parameters 


Parameter 

Value 

Unit 

Ao 

10 

V 

A, 

2 

V 

A2 

2 

V 

fi 

0.9 

GHz 

h 

1 

GHz 

R 

1,000 

Ohm 

k 

1.00E-16 

Amp 

Table  1 :  Circuit  And  Input  Parameters. 

The  Arithmetic  Operator  Method  requires  the  input  signal  to  be  a  spectral  vector  of  the  form 
X  =  [Xo  Xu  X|i  Xu  Xji . .  .Xk,  Xkif .  where 

•  Xo  is  the  dc  term, 

•  k  is  number  of  frequencies  in  the  BIPD  table,  and 

•  Xicr,Xki  represents  the  complex  number  pair  (real  part  and  imaginary)  describing  the  amplitude 
and  phase  of  the  kth  frequency  component. 

To  reach  this  form,  rewrite  (1.41)  as 


A  A 

Vi=Aq[\-^x],  where  x  =  —  sin(2;ry;0  +  — sin(2;r/20  (1.42) 

^0  ”^0 

Now,  for  the  input  signal  in  (1.42),  the  AOM  input  is  the  frequency  domain  spectral  vector  ofx,  X,  and 
consists  of  the  following  components: 


Xo=0 
X„  =  0 
Xii  =  -A./  Ao 
X2r  =  0; 

X2i  =  -A2/  Ao 

As  will  be  explained  further  in  just  a  moment,  there  are  10  frequencies  in  the  BIPD  table,  so  the  spectral 
vector  X  must  contain  2*10+1  ==  21  terms.  Thus,  to  form  the  spectral  vector  X  simply  pad  the  input  vector 
with  zeros  to  obtain  the  correct  size.  The  spectral  vector  for  this  simulation  is: 

X  =  [0,0, -(A,/  Ao),0,-(A2/  Ao), 0,0, 0,0, 0,0,0, 0,0, 0,0, 0,0, 0,0, of  (1 .43) 

Now  we  must  create  the  Spectrum  Transform  Matrix  T* .  To  do  this,  we  need  to  choose  a  set  of  BIPD 
frequencies.  For  this  example,  we  are  interested  in  all  of  the  frequencies  in  the  frequency  domain 
surrounding  the  fundamentals  and  its  second  harmonics.  We  could  choose  to  have  more  frequencies,  but 
since  we  can  sometimes  put  second  harmonics  to  good  use  [12],  while  third  harmonic  content  is  of  little  use 
and  can  be  filtered,  we  will  choose  the  second  harmonic  content  limit.  This  gives  us  the  following  BIPD 
table: 
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kx 

BIPD  Table 

ni 

Pa 

1 

1 

0 

2 

0 

1 

3 

-1 

1 

4 

2 

-1 

5 

-1 

2 

6 

2 

0 

7 

0 

2 

8 

3 

-1 

9 

1 

1 

10 

-1 

3 

Table  2:  BIPD  Table. 


Notice  in  Table  2  that  we  are  not  ignoring  third  order  intermodulation  products,  as  the  indices  for 
=  4and  =  5  show.  Furthermore,  there  will  be  some  fourth  order  intermodulation  products  that 

appear  in  the  range  of  the  second  hamionics,  so  we’ve  added  the  frequency  indices  =  8  and  =  10  in 
order  to  take  account  of  this.  To  gain  some  visual  insight  into  the  set  of  frequencies  we  are  considering, 
the  Table  2  above  is  plotted  graphically  in  Figure  2  below  using  our  selected  values  for  the  non- 
commensurate  tones,  =  0,9  GHz  and  /j  =  1 .0  GHz  .  (Note:  The  plot  shown  in  Figure  2  accurately 
shows  the  locations  of  the  BIPD  frequencies,  but  does  not  give  an  accurate  depiction  of  the  output 
amplitudes  we  can  expect  from  the  logarithmic  amplifier.  This  is  because  the  plot  was  done  assuming  an 

arbitrary  transfer  function  of  the  form  y(r)  =  x(t)  +  (/)  +  a:’(0  + (0  with 

x(t)  =  Ibos(27l  fjt)  +  2!bos(2;r fy)  used  as  the  input,  and  where  some  of  the  fourth  order  mixing 

components  have  been  intentionally  dropped.) 
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Frequency  (BIRD  Indices  indicated) 


Figure  2:  BIPD  Locations. 

Comparing  the  table  above  with  that  shown  earlier  in  the  2.6.1  Basic  Intermodulation  Product  Description 
Table  section,  we  do  not  carry  DC  as  a  BIPD  in  the  Table  2.  We  did  not  overlook  DC  in  our  AOM 
modeling,  but  we  have  found  that  it  is  more  convenient  to  take  DC  mixing  terms  into  account  at  the  time 
the  Spectrum  Mapping  Table  is  constructed  when  using  computer-aided  methods,  and  this  is  what  we  have 
done.  This  is  because  the  DC  terms  created  through  mixing  of  BIPD  entries  are  essentially  composed  of 
fundamental  inputs  (designated  hy  and  =  2  in  the  table  above)  mixing  with  themselves  through 

subtractive  intermodulation.  Using  the  systematic  methods  described  previously,  the  Spectrum  Mapping 
Table  is  created,  and  from  this  the  Spectrum  Transform  Matrix.  For  the  example  here  with  10  non-zero 
BIPD  terms,  the  Spectrum  Mapping  Table  has  115  entries.  With  the  Spectral  Vector  and  the  Spectrum 
Transform  Matrix  now  defined,  we,  perform  the  series  expansion  of  (1.36)  in  the  frequency  domain  for  the 
output  spectral  vector  Y : 


Y  =  X-lT,X+iT,(T,X)-lT,(T,(T,X))  +  ...  (1.44) 

For  this  simulation,  the  infinite  series  expansion  as  described  in  equation  (1 .36)  was  truncated  at  20  terms. 
Figure  3  below  shows  the  accuracy  of  the  20  terms  approximation  over  the  interval  (—  1  <  X  <  1)  .  Note 
that  with  the  magnitude  x  of  the  input  having  a  maximum  value  of  0.4,  we  can  expect  veiy  little  error  in  the 
results  attributable  to  the  decision  to  truncate  the  series  expansion  at  20  terms.  But  since  the  number  of 
terms  has  a  direct  impact  on  simulation  time,  minimizing  the  number  of  terms  given  an  accuracy  target 
should  always  be  a  goal.  We  established  our  accuracy  target  through  graphical  inspection  of  Figure  2,  but 
it  is  feasible  to  dynamically  adapt  the  number  of  temis  depending  on  the  iteration  error.  This  is  a 
convenient  place  to  optimize  a  simulation  beyond  what  is  presented  here. 

Transforming  Y  back  to  the  time  domain,  and  adding  the  constant  terms  of  (1 .40),  the  AOM  simulation 
results  compare  favorably  with  Spice  output  as  shown  in  Figure  4. 
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The  results  indicate  good  agreement  with  the  Spice  output.  The  small  discrepancies  that  one  sees  can  be 
attributed  to  two  factors:  First,  the  Spice  model  used  a  gain  of  100  for  the  operational  amplifier  rather  than 
an  infinite  value,  so  the  Spice  output  itself  may  have  some  error  compared  to  a  completely  ideal  model. 
Running  Spice  with  the  op  amp  gain  as  a  parametric  value  would  make  for  an  interesting  area  of  further 
research.  Second,  since  we  have  demonstrated  that  a  20  term  polynomial  expansion  for  the  natural 
logarithm  function  yields  an  accurate  approximation  given  the  magnitude  limits  on  x,  the  error  is  likely 
attributable  to  deviations  in  the  Spectrum  Transform  Matrix.  We  get  deviations  in  the  Spectrum  Transform 
Matrix  through  our  choice  of  BIRD  frequencies.  It  is  likely  that  increasing  the  number  of  BIRD  frequencies 
will  lead  to  more  accurate  results.  Referring  back  to  Table  2,  we  chose  to  limit  our  fourth  order 
intemiodulation  terms  to  only  two  -  those  for  =  Sand  =  10 .  We  will,  however,  also  have  fourth 

order  intermodulation  products  at  DC  due  to  subtractive  mixing  of  the  second  harmonic  frequencies,  these 
were  not  taken  into  account,  and  it  is  possible  that  their  inclusion  could  increase  the  accuracy  of  the  results. 


4.  Conclusion 

We  have  explained  the  Arithmetic  Operator  Method  (AOM)  mathematics  in  some  detail  and  shown  the 
relationship  between  the  mathematics  and  the  algorithms  used  to  implement  the  method  in  the  paper  by 
Chang  and  Steer  [2].  We  have  also  presented  a  case  study  of  the  application  of  AOM  for  a  logarithmic 
amplifier  and  demonstrated  that  the  results  of  using  the  method  compare  favorably  to  those  from  a  Spice 
reference  simulation  model.  It  is  our  view  that  these  favorable  results  confirm  the  utility  of  the  method  for 
behavioral  modeling  of  devices,  and  it  is  our  hope  that  our  explanations  and  case  study  will  help  others 
apply  AOM. 
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